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Abstract 

Traditional numerical techniques for solving time-dependent partial- 
differential-equation (PDE) initial-value problems (IVPs) store a trun- 
cated representation of the function values and some number of their time 
derivatives at each time step. Although redundant in the dx ^ limit, 
what if spatial derivatives were also stored? This paper presents an iter- 
ated, multipoint differential transform method (IMDTM) for numerically 
evolving PDE IVPs. Using this scheme, it is demonstrated that stored 
spatial derivatives can be propagated in an efficient and self-consistent 
manner; and can effectively contribute to the evolution procedure in a 
way which can confer several advantages, including aiding solution ver- 
ification. Lastly, in order to efficiently implement the IMDTM scheme, 
a generalized finite-difference stencil formula is derived which can take 
advantage of multiple higher-order spatial derivatives when computing 
even-higher-order derivatives. As is demonstrated, the performance of 
these techniques compares favorably to other explicit evolution schemes 
in terms of speed, memory footprint and accuracy. 

1 Introduction 

In traditional methods for solving time-dependent partial-differential-equation 
(PDE) initial- value problems (IVPs) , a truncated representation of the function 
values and some number of their time derivatives, on some number of spatial 
slices, are stored on a discrete grid. As is well known, the grid may consist 
of regularly-spaced points or may be an irregular point collection, and there 
are many choices for the basis used to represent the function values at the 
grid points. Nevertheless, evolving the PDE inevitably requires reconstructing 
spatial derivative values, and because these values are not stored explicitly, 
they must be reconstructed from the truncated representation of the function 
valuesj^ This paper reports on an investigation into the construction of a scheme 

^For spectral methods that use a Fourier basis, this reconstruction is trivial. Unfortunately, 
spectral methods are often impractical for problems with nonlinearities, irregular geometries 
or nonperiodic boundary conditions. 



which, unlike traditional schemes, stores some number of spatial derivatives on 
the grid along with the function values. Although redundant in the dx — >■ limit, 
at least in the context of the scheme presented here, storing spatial derivative 
values will be shown to have practical advantagetj^ To construct a scheme that 
can make use of, and evolve, stored spatial derivatives, it is natural to consider 
methods based on generating local power-series expansions of the solution. 

Constructing power-series solutions to ordinary differential equations is a 
basic technique included in most introductory texts on mathematical methods 
or differential-equations (for example, see [4]). Algorithmically constructing 
power-series solutions to ordinary differential equations has been extensively 
studied (see [3], and references therein), but has remained a little-known tech- 
nique. Although sometimes used for the numerical evaluation of special func- 
tions [D, the construction of power-series solutions has been generally thought 
of as an analytic tool and not as the basis for numerical algorithms. This 
is changing, and algorithmically constructing power-series solutions to ordi- 
nary differential equations is gaining in popularity. This is now often called 
the Differential Transform(ation) Method (DTM) [S], but is also known as the 
Parker-Sochacki method [20], or generically as the Picard-iteration method or 
the Taylor-series method. In this paper, I adopt the notations and terminol- 
ogy from the DTM literature. The conventions which define the DTM can 
be straightforwardly extended to multivariate power series for application to 
PDEs, but applying the resulting definitions to solve PDEs numerically requires 
dealing with several complexitie^ This paper details a practical method for 
numerically applying the DTM to time-dependent PDEs: an iterated, multi- 
point DTM (IMDTM), and shows that the resulting evolution algorithm pos- 
sesses several highly-desirable properties. Notably, by storing a sufficiently-large 
set of spatial derivatives, the resulting system can be configured to possess a 
directly-calculable self-consistency constraint, the violation of which measures 
the quality of the solution. 

The remainder of this paper is organized as follows: Section [2] reviews how 
the evolution equations for higher-order derivatives can be derived. Section|3]re- 
views the multivariate DTM formalism. The IMDTM algorithm is presented in 
Section [4] The efficient implementation of IMDTM requires an efficient polyno- 
mial interpolation scheme capable of using higher-order derivative information. 
Such a scheme is also presented in Section |4] The paper then concludes, and a 
set of appendices follow. 

There has been hmited investigation of using stored spatial derivatives in the context of 
finite-volume schemes for evolving systems of conservation laws. In the finite-volume litera- 
ture, these are known as multi-moment schemes (see [ISJ and references therein). 

•^In the context of finite- volume schemes, using local (multivariate) power-series expansions 
to enhance the accuracy of Riemann solvers forms the basis of the well-known Arbitrary 
accuracy DERivative (ADER) Riemann technique 1221 . Unlike ADER schemes, the technique 
presented in this paper can be applied to any kind of time-dependent PDE. 
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2 Evolution Equations for Higher- Order Deriva- 
tives 



Before discussing the conventions which define the DTM, consider the following 
two examples which explicitly demonstrate how the higher-order derivatives of 
a PDE's solution obey PDEs derivable from the PDEs for the solution function. 
First, consider the 2-D wave equation: 

from which the evolution equation for the higher-order spatial derivatives can 
be derived by applying the operator 

t =0 (2) 

Because the wave equation is linear, all of the derivatives obey structurally- 
identical equations. For nonlinear equations, this is not the case. For example, 
consider the modified KdV equation: 

for which the evolution equation for the higher-order spatial derivatives are: 

Wtd^^^^\k) 1 ^ \j) d^ dx^-i j ax"-'=+i '^dx^d^~ ■ ^' 

To avoid a "combinatorial explosion" in the number of terms, I have used (mul- 
tivariate) versions of the Leibniz rule (the product rule) and the Faa di Bruno 
formula (the generalized chain rule) [10]. The multivariate DTM encapsulates 
these formulas in easy-to-apply recursion relations, thus avoiding the direct 
use of combinatorial calculus formulas as was done in Equation |4j Because the 
higher-order coefficients depend on the lower-order coefficients, the higher-order 
coefficients will tend to be slower to compute, and in some cases less numerically 
stable, than the lower-order coefficients. 

It is possible to evolve all of the higher-order coefficients using a traditional 
PDF evolution scheme by discretizing the associated evolution equations, which 
will, generally, increase the computational time required by at least a factor of 
n for evolving n spatial derivatives. The IMDTM scheme allows multiple coeffi- 
cients to be evolved together, significantly reducing the computational overhead. 



3 Multivariate Differential Transform Method 

In many cases, the application of the classic integral transform methods (i.e. 
Laplace and Fourier) can be reduced to use of a table of substitutions . Simi- 
larly, a differential equation can be transformed into a recursion relation for the 
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coefficients of its power series solution(s) using a table of substitutions [5], and 
recent literature refers to this as The Differential Transform Method [5]. The 
conventions which define the DTM can be extended to the multivariate case in 
the obvious way specifically: 



W{ki,k2,...,kn) 



1 



so the inverse transform is: 



Qki+k2-{ hfc„ 



0,0,. ..,0 

(5) 



w{wi,W2, ■ ■ .,W„) 



EE 

ki=0 k2=0 



E ^(^1 



)n 



(6) 



and the basic operations are given in Table [T] By combining Equations [5] and [6] 
it is clear that the multivariate DTM is nothing more than the construction of a 
multivariate Taylor expansion. The efficient evaluation of nonlinear terms is an 
important practical concern. Fortunately, efficiently computing nonlinear func- 
tions of power series is a well-researched problem by the creators of automatic- 
differentiation software, and is generally done either by directly evaluating a set 
of multivariate recurrence relations |23j or by combining different evaluations 
of univariate recurrence relations using different directional projections of the 
original series [2] (TH] . Direct evaluations of the multivariate recurrence relations 
will be used here, and the DTM translations for additional (nonlinear) functions 
are provided in Appendix [B] 



Original Function 


Transformed Function 


w{x) — y{x) ± z{x) 

w{x) — Xy{x), A a constant 

w{x) = y{x)z{x) 
wlx) = 1 • • • x;^" 


W{k) = Y{k) ± Z{k) 
W{k) = \Y{k) 

w{k) = ('^'+r^'y(fc + r) 

W{k) ^Y!ll^,---Y!l:^,Y{l)Z{k-l) 

W{k) ~ nr=i ^ki,mi, S is the Kronecker delta 



Table 1: Basic operations under Multidimensional DTM 

Note that in the example 1 -I- 1-dimensionaQ systems discussed in later sec- 
tions, the temporal coordinate is called t and the spatial coordinate is called x. 
Correspondingly, k is the DTM-coefficient index for temporal derivatives, and 
h is the DTM-coefficient index for spatial derivatives. Also, an overdot will be 
used to indicate a time derivative. 

*l-|-l-dimensional means that the system has one spatial dimension and one time dimen- 
sion. 
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4 Iterated, Multipoint DTM 



Although the multivariate DTM has been applied as a convenient way to analyt- 
ically construct a power-series solution to PDEs [3] [14] , a novel contribution of 
this paper is a method for using the multivariate DTM as a numerical evolution 
scheme which can handle an arbitrarily-large Cauchy surface. When solving 
ODEs, the initial conditions can be exactly represented using a finite set of val- 
ues. For PDEs, however, the initial conditions are themselves (differentiable) 
functions and they cannot be represented, in general, using only a finite set of 
numbers. As expected, applying the multivariate DTM to a PDE results in a 
recurrence relation which allows higher-order derivative coefhcients in one vari- 
able to be calculated in terms of the complete set of derivatives with respect to 
the remaining variables (i.e. on the Cauchy surface). For example, consider the 
2-D wave equation: 

for which the DTM gives the recurrence: 

m'.)= "'y,'i\;" m-2,;. + 2). (8) 

Given the initial tower of derivative coefficients, F{0, h) and F{1, h) (the spatial 
derivatives of / and /, respectively), any F{k, h) can be calculated. 

Even if the radius of convergence of a power series is infinite, when working 
at finite precision, the useful range may be limited, even as more terms are used. 
As a result, the IMDTM uses power-series expansions around multiple points on 
the Cauchy surface to precisely represent the initial conditions. In other words, 
at each point on the spatial hypersurface, a tower of DTM coefficients are stored, 
{F(fc, h)t^x}- For a PDE which is first-order in time, only the fc = coefficients 
are stored; for an equation which is second-order in time, only the k E (0, 1) 
coefficients are stored; and so on for higher-order PDEs. In combination with 
the recursion relation derived from the PDE, the power-series expansion of the 
solution at each spatial point can be reconstructed from the stored coefficients. 
This provides the IMDTM scheme with a powerful feature: it can yield an inter- 
nal self-consistency condition. Specifically, the power series around some point 
must, by self-consistency, provide the values of the function at the neighboring 
points. As the system is evolved forward in time (via iteration), this condition 
will break down, and this breakdown can be used to measure the quality of the 
evolved solution. 

As you might expect from the definition of the Taylor expansion, a time step 
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from t = to t = dt is calculated as: 
F{0,0)t=dt,x= f{x,dt) = 

F{0,l)t^dt,.^ ^f{x,dt)^ 



Y.F{k,0)t=o,.dt'' (9) 

oo 

J2{k + l)F{k + l,0)t=o^,dt'' (10) 

fc=0 

oo 

Y,F{k,l)t=o,.dt'' (11) 



fc=l 



where, again, the overdot indicates a time derivative. 

At any time, i, the system is specified by an infinite set of coefficients, thus 
the description must be truncated for numerical implementation. This trunca- 
tion immediately generates an impediment to the construction of an iterative 
scheme: at whatever order the truncation is performed, even-higher-order terms 
are necessary to compute the evolution of the highest-order terms being stored. 
If these missing coefficients are simply taken to be zero, then the scheme will 
quickly destabilize. Qualitatively, if n is the highest order stored, and all terms 
of an order greater than n are taken to be zero, then the term of order n will 
remain constant, the term of order n — 1 will grow linearly with time, the term 
of order n~2 will grow quadratically with time, and so on. For the scheme to be 
practical, a method for computing these missing coefficients must be provided. 



4.1 Interpolating Polynomials 

In many traditional grid-based numerical evolution schemes, the derivatives are 
approximated by a finite-difference computation. Conceptually, this approxi- 
mation uses some number of neighboring values to construct an interpolating 
polynomial, and uses the derivative of that polynomial (at the central point) 
to estimate the derivative of the function represented on the grid [24J. Since 
the IMDTM stores truncated spatial power series on the grid, instead of just 
the function values, it is possible to use multiple coefficients per point (the 
lower-order derivative values) in order to reconstruct the needed higher-order 
derivative values. However, there are several complications which need to be 
discussed. 

First, because the coefficients are stored only to some fixed precision, not 
all of the neighboring-point power series coefficients should be used for the 
derivative-reconstruction procedure. To understand this, assume that the power 
series at a point can be used to calculate the function value at any neighboring 
point to some precision, and that the stored function value at any neighboring 
point is accurate to that same precision. This means that at least the neighbor- 
ing zeroth-order components have no significant information to contribute to 
the higher-order derivative interpolation. To put it another way, since it takes 
n orders (i.e. n terms of the power series) to converge to the desired precision. 
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that limited precision is equivalent to introducing an error term of order n + 1. 
Adding numbers with an error of order n + 1 to calculate an order- (n + 1) coef- 
ficient would be impossible. Even restricting the interpolation procedure to use 
only the last few orders, it is still, in general, not possible to reconstruct the 
missing derivatives to the same precision as the inputs. Luckily, by making dt 
small enough, reconstructing the missing higher-order coefficients to the same 
precision as the inputs is unnecessary. 

The naive method for constructing a polynomial interpolant using higher- 
order derivative information - by solving the implied matrix equation - is not 
practical for higher-dimensional systems. Fortunately, it is possible to write 
down an explicit solution to the general interpolation problem, and use that 
solution to calculate the required coefficients. The method used here, detailed 
in Appendix |Aj uses the multipoint Taylor expansion as defined by Lopez and 
Temme |18j to explicitly construct the interpolating polynomial. Although the 
method seems relatively complicated, a number of the needed factors depend 
only on the geometry, specifically the relative distances between neighboring 
grid points, and can be cached once calculated (or, precalculated, if the prob- 
lem geometry is static). After all of the geometry-dependent factors have been 
calculated, a linear combination of the input coefficients from the various neigh- 
boring points yields the missing higher-order coefficients. This method can be 
thought of as a generalization of a finite-difference stencil that can make use of 
higher-order derivative information. 

An IMDTM scheme where many higher-order coefficients are stored on the 
grid and interpolation is used only for the highest-order coefRcients is numer- 
ically unstable. This arises for the same reason that the missing higher-order 
coefficients cannot be set to zero: a constant error in the order-n coefficient will 
lead to an error in the order- (n — 1) coefficient which grows linearly with time, 
and so on. In practice, only a couple of orders can be evolved together. You can 
store only a couple of orders on the grid and use interpolation to generate the 
coefficients necessary for forward evolution. That can be stable, but lacks the 
self-consistency constraint. A stable scheme with a self-consistency constraint 
can be constructed by "stacking" a number of such lower-order schemes on top 
of one another. For example, interpolate using orders zero and one to evolve 
orders zero and one, interpolate using orders two and three to evolve orders two 
and three, and so on. This is illustrated in Figure [l] The following examples 
demonstrate this technique. 

4.2 Example: The Wave Equation 

The salient features of the IMDTM scheme are highlighted with a simple exam- 
ple: the l-|-l-dimensional linear wave equation. This example, like the nonlinear 
equation example which follows, have been chosen because they have known 
periodic, non-singular, analytic solutions. The implementation uses periodic 
boundary conditions so that the interpretation of the results is not complicated 
by details of more-complicated boundary-condition handling. Since each grid 
point stores an entire tower of derivatives, the boundary points would need to do 
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this as well. How best to use this freedom for other kinds of boundary conditions 
will be the subject of future research. 

For this demonstration, the IMDTM implementation used a grid of = 16 
points with the coefficients of the first 14 spatial orders, including the constant 
term, stored at each point for both the field value and its first time deriva- 
tive (i.e. F(0,0) through F(0, 13) and F{1,0) through -F(l, 13) are stored 
at each point). The physical size of the grid is L = 18, so dx = 1.125 and 
dt = 1 is used for the temporal evolution. In this case, the values of dx and 
dt are exactly-representable floating-point numbers, although the accuracy is 
not significantly degraded if this is not the case. The initial conditions were 
^0(2;) = cos(27rx/L), 0g(a;) — 0. Unlike traditional explicit temporally-higher- 
order algorithms (e.g. RK4), IMDTM does not require additional grid copies to 
hold intermediate grid state^ Nevertheless, the performance of the IMDTM 
implementation is compared to an RK4-based scheme. For the particular sys- 
tem chosen here, a second-order stencil (with dx ~ 0.08) is accurate to approx- 
imately 4 significant digits, and an eighth-order spatial stencil is accurate to 
over 13 digits. The dt values for the RK4 scheme were chosen to match the 
single-core running time to that of the IMDTM implementatiorj^ 

Next, consider the most-basic stable use of the IMDTM scheme: only the 
function values (the zeroth-order coefficients) are stored on the grid. The erroi[^ 
of this method in solving the wave equation is illustrated in Figure [2] In such a 
configuration, the scheme becomes a kind of fancy finite-differencing algorithm. 
As expected, using a larger interpolation radius leads to a more-accurate evolu- 
tion; this is exactly equivalent to using a higher-order finite-differencing stencil. 

For the case where only the function values are stored, and using a unit 
radius for interpolation (meaning each point and its two immediate neighbors 
contribute to the interpolation), it is practical to write down the explicit time- 



stepping rule. Expanding Equations 26 and |30| (using Equations 31 and 32 as 
noted)j^ and using the recurrence relations derived from the wave equation's 
PDF (given in Equation |8| with Equations |9] and 10 yields the following time- 

^The method of (spectral) deferred corrections can also be used to produce higher-order 
temporal evolution for general time-dependent PDEs without storing several intermediate 
grid copies [?]■ As mentioned previously, for systems of conservation laws, ADER schemes 
can generate higher-order temporal evolution without storing intermediate grid copies. 

^AU code was compiled using GNU gH — h version 4.4.4 given the flags: -03 -march=native 
-msse -mfpmath=sse -DNDEBUG. 

'^The average error at each time slice is computed by taking the mean of the base-ten 
logarithm of the relative error at each point. 

*The equations in Appendix |a] use a to represent the multiset of spatial derivative orders; 
to apply those formulas to the l-fl-dimensional case here, note that a always has unit cardi- 
nality, and its only integer element runs from zero to infinity. In practice, the sum must be 
truncated at some finite order (determined by the maximum possible order of the interpolat- 
ing polynomial, which is determined by the number of neighbors and the number of stored 
spatial derivative orders used for interpolation). 
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stepping rule: 



f{t + dt,x) = ^^{-18{dt^-2dx^)f{t,x)+ 

dt (9dtf{t, x~dx)+ 9dtf{t, x + dx) - 2dt^f{t, x) + 
mdx^j{t, x) + dt^f{t, x-dx)+ dt^f{t, x + dx)j ^ 

f(t + dt,x) = — ^(-12dtf(t,x) + 6dtf(t,x-dx) + 6dtf(t,x + dx)- 
12dx^ 

2dfj{t, x) + 12dx^f{t, x) + dt'^fit, x - dx) + dfj{t, x + 

where the overdot indicates a time derivative. Note that because the interpo- 
lating polynomials are fitting only one value from each of three points, they are 
only quadratic functions, and so can only estimate the first and second spatial 
derivatives at each point. So, at each point, {F(0, 0), 0)} are given and 
{F(0, l)^^i^(0, 2),F(1,1),F(1,2)} are estimated using the interpolating poly- 
nomialq^ From the recurrence relation given by Equation Isl the values of 



{i^(2, 0), i^(3, 0)} can be calculated. The sums in Equations P and 10 can be 
truncated accordingly. Because the PDE is linear, these same equations are 
obeyed for each spatial derivative order. The point of this paper, which is 
that there can be a benefit to storing and evolving spatial derivative values, is 
demonstrated next. Figure [3] shows the case where two orders are stored on the 
grid and both of those orders are used for the interpolation procedure. Even 
with a smaller interpolation radius, the result is much more accurate than in 
the preceding case. As expected, using more coefficients from closer grid points 
yields a more accurate result. 

The IMDTM implementation yields very-high precision evolution while tak- 
ing large time stepip°| This could be advantageous for parallel implementa- 
tions on non-shared-memory clusters, because each time step involves slow 
synchronizing communication. As for the construction of a scheme with the 
aforementioned neighbor-to- neighbor self-consistency constraint, the suggested 
mechanism of "stacking" evolution schemes for the coefficients in order pairs 
is demonstrated in Figure |4j As can be seen by comparing the neighbor-to- 
neighbor constraint-violation error and the analytic error, the solution is more 
self-consistent than it is accurate. Comparing the absolute errors instead of the 
relative errors does not change this conclusion. 



^To be clear, two interpolating polynomials are constructed: One to estimate the spatial 
derivatives of f{t,x) and one to estimate the spatial derivatives of f{x,t). 

^"^This linear case does not form a fair basis for benchmark-like comparison with other 
numerical evolution schemes, such as an RK4 implementation, in general, because evaluating 
the recursion relations necessary for the computation of the nonlinear terms can add significant 
expense, whereas the same is not true for schemes which work only with the function values. 



9 



4.3 Example: A Strongly-Nonlinear Equation 

While solving the linear wave equation allows the demonstration of several fea- 
tures of the IMDTM scheme, it is not restricted to linear problems. As an ex- 
ample, the IMDTM scheme can be used to evolve the strongly-nonlinear PDE 
in Equation [sj This equation has a non-singular, periodic solution |16j : 

-^^"+ 2±cos(2ax-8«3t) - ^^^^ 

where a is a free real parameter. Applying the multivariate DTM (see Table [I]) 
to Equation [3] yields: 

fc h 

H{k,h) = ^^F{k-m,n)F{m,h-n) (15) 

m— n— 
k h 

G{k,h) = ^ ^(/i- n-t- l)i7(fc- m,n)F(m,/i- n-Kl) (16) 

m— n— 

F(k,h) = -I {G(k - 1, h) + (h + 3)(h + 2)(h +l)F{k- I, h + 3)) (17) 
k 

where H{k, h) corresponds to the p factor, and G{k, h) corresponds to the p g| 
term. Caching the computed values of i?(fc, h) is essential to an efficient imple- 
mentation. The stable evolution of Equation [M] as the solution of Equation |3] 
requires a smaller dx and a smaller dt compared to those used for the linear wave 
equation. For this demonstration, the IMDTM implementation used a grid of 
N — 78 points. The physical size of the grid is i = 43.875, so dx = 0.56250 
and dt — 0.001 is used for the temporal evolution. The initial conditions were 
(j)o{x) — -~2\f2a -I- 2+m&i^ax) ' ° ^ tt/L. As should be expected given the range 
of the sums in the recurrence relations, unlike for linear PDEs, computing the 
higher-order terms is significantly more expensive than computing the lower- 
order terms. Figure [5] shows how the IMDTM scheme, when only two orders 
are stored on the grid, compares to an RK4 implementation performing the same 
amount of computational work and using the same number of grid degrees-of- 
freedom (not counting the auxiliary grid copies used for the intermediate RK 
steps). However, reducing the RK4 time step to dt = 0.001, identical to that 
used for the IMDTM scheme, does not significantly degrade the RK4 accuracy. 
While the IMDTM scheme is less stable than the RK4 scheme, it is far more 
accurate for tens of thousands of time steps. Although the IMDTM scheme does 
need extra buffers to efficiently evaluate the nonlinear terms, for this equation 
it still uses significantly less memory than the RK4 implementation. Although 
it might seem "more fair" to compare to a RK16 scheme, the memory use of an 
RK16 implementation would be almost 10 times that of the IMDTM scheme, 
and would likely be unsuitable for practical, large-scale implementation. 

Pair-wise order evolution (as illustrated in Figure [ij was used in order to 
generate a scheme with a neighbor-to- neighbor self-consistency constraint, and 
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the results were similar to those for the hnear case (shown in Figure |6]). The 
growth of the analytic error over the neighbor-to-neighbor constraint-violation 
error is exacerbated when computing at higher precision. SpecificaUy, when 
using the double-double type from Bailey's qd package |T3j, which has twice the 
precision of an IEEE double, after some hundreds of time steps, the analytic 
error becomes greater than the neighbor-to-neighbor constraint-violation error. 
This is demonstrated in Figure [Tj 

5 Future Work 

The work presented here leaves a lot of room for future research. The following 
are some of the items which deserve further investigation: 

• A procedure for performing a stability analysis of an IMDTM scheme 
should be developed. 

• The IMDTM should apply naturally to problems with nontrivial boundary 
conditions, so long as the multiple power series expansion can be generated 
in a self-consistent manner. The exact conditions under which this is 
possible should be investigated. 

• DTM can be applied to boundary-value problems by leaving some of the 
lower-order coefhcients free and then solving for them in a way consis- 
tent with the boundary constraints after expressions for the higher-order 
coefficients have been obtained [11] [12]. Similarly, it may be possible 
to extend IMDTM to handle multivariate boundary-value problems and 
evolution problems with constraints. 

• Alternate methods for computing the power series at time t + 5t from 
those at time t which are more stable than the method presented here 
might exist. For example, it has been demonstrated that using Fade 
approximants can increase the stability of an iterated DTM scheme [21] [9] . 
It may also be possible to formulate the IMDTM as an implicit scheme. 

• A reusable framework for developing IMDTM codes is under development, 
and its interface definitions and features will be presented in a future 
publication. 

6 Conclusion 

In this paper, the IMDTM PDE evolution scheme for initial- value problems has 
been introduced. Methods to stabilize the evolution and efficiently use polyno- 
mial interpolation to enable iterative forward propagation have been detailed. 
It is an explicit evolution scheme, and it tends to be able to take larger time 
steps compared to other explicit schemes (e.g. RK4) while performing a simi- 
lar amount of computational work. For many equations, the IMDTM scheme 
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will also use less memory than the corresponding RK-like scheme. The larger 
time steps should allow for substantially-reduced communication on non-shared- 
memory machines. Because the algorithm is local and explicit, it should scale 
strongly to the largest conceivable siipcrcomputers. Fiirthcrmore. IMDTM can 
naturally be used to construct a verifiably-self-consistent code in the sense that 
it has an internal self-consistency constraint which can be used to measure the 
quality of the solution. Not only does this often allow the user to know when 
the solution can no longer be trusted, it can allow the IMDTM solution to be 
used as a baseline for verifying other codes. More generally, this paper demon- 
strates that efficient numerical schemes can be constructed which store and 
evolve higher-order spatial derivative information, and that doing so can confer 
significant practical advantages. 
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A Interpolating Polynomial 



The interpolating polynomial can be derived from the multipoint Taylor expan- 
sion as defined by Lopez and Temme [18]. Let {x[^\x\'^\ . . . be the set 
of m unique coordinate values in dimension i € {1, • • • ,D}. In one dimension, 
a function f(x) has the multipoint expansion: 



CO I m TT™ 



fix) - V I V nLi,fc^,(a^ 



\[{x-x^^'^T 



k=l 



(18) 



where 



E 



1 



Jn — 1 



1 d" 

n! dx" 



(n - 1)! 



fc=l,fc/j 

which can be written as: 

_ " 1 d"^(M3) 

~ (71 - (fc ^ j))! dx"~(k^3) 



(19) 



fe 



This can be generalized in the usual way to a multidimensional function by 
recursive expansion, yielding: 



n T— rm / ^, 

D rn 



(21) 



j=i k=i 

Expanding the final, non-constant binomial term gives: 



a \ ^e{l,...,m}0 



7e{0,...,m}x{0,l}'"^M7|=rr 



70 T-rni ('„^('^')\7fc 



D 



xnn E 

i=\ k=\ \p=0 



(-xf ))--^'< . (22) 



Note that the product of two sequences of length N is their convolution: 



AT 



N 



2N 



E E = E E ^" 



(23) 



i=0 \j=0 
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and the same applies for three sequences: 



N 



N 



N 



'2N / n 



^ N 



\n=0 



Kn=0 



\n=0 



Vi=0 \i=0 / / \n=0 / 



from which the pattern can be seen. 



/(^) = E ( E n^^l ET-g{ 0,...,m}x{0,l}">-M7| 



70 1-rm C-T^ M 



i 

D fmOLi p 91 



n E<EE-E r K-i"^) 



i=l \p=0 91=0 92=0 9m=0 



„(™)^Q,-9,r 



V^m-l - 9m/ \P- qij 



i-(p-9l) 



The contribution to the DTM coefficient F{h) is then easily derived from the 
coefficient of f{x): 



/(a;)-E E T-rD ,-rm / 



_D m a; p 

xHE^E 



E n (- 

i=l p=0 j=0 Y7e{oa}™-i,|7|=m-l-(p-j) fc=l,fc5^/3i 

xfEE---E f:0(--^^)"'"" 



^91=0 92=0 9m=0 



OLi 

1m -1 — In 



(_2,(™-l))a>-(9™-l-9™) . . . ^ j (^^x^^^-'jai-U-qi) 



(26) 



and 



E 



7e{l,...,m}° 



1 

a! 9a;" 



.(ft) 



(27) 

where a — a — {ji ^ /3i}. The product rule for a general partial derivative 
is HOl: 



^fciH hfc. 



i;:M = E---E 



(28) 
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So we can expand the expression for 



7e{i, 



a\ ^ — ' \6 J dx^ 



and the factor for each dimension can be separated: 



0L\ '■ — ' \ (5 / dx^ 

7e{l,...,m}-° (5,Vii5i<Qi ^ ' 



dx', 



i-<5. 



Because ^ In / = i£ we can write ^ = In/, and so: 



(29) 



(30) 



(ft) 



(fe) 



The second factor in Equation 30 can be evahiated recursively using Equation 31 
and the product rule: 



dx' 



-uv 



E 



fc=0 



n\ d" ''u d'^v 



(32) 



J dx" ^ dx'^ 

B Additional Multivariate DTM Recurrences 

The following table shows how to apply the Multidimensional DTM to some 
common nonlinear terms. Once the formulas for the basic arithmetic operations 
are known, deriving the formulas for functions that are solutions of an ODE is 
straightforward [8] [17] (just apply the DTM transformation rules to the defining 
ODE). A large table of these recurrences can be found in an appendix of the 
Fast Forward Automated Differentiation Library (FFADLib) User Manual ^5] . 



Original Function 


Transformed Function 


w(x) ~ y{x)/z{x) 
w{x) = \/y{x) 


Wik) = ^ 




i)] 
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'w{x) = exp(y(a;)) 
w{x) = \n{y{x)) 

w{x) = y{xy 

w{x) = sm{y{x)) 
w{x) = cos{y{x)) 



w^(fc) = ^E'U---Ef:Jo'---EL=o 

{ka-la)W{l)Y{k-l) 



W{k) 



Y(0) 



ia=0 ' 



E, 



ika-la)Yil)W{k-l) 



W{k) 



Y(0) 



_.T4^(o)y(fc) + ^E?:U---Et"=o 

(fc„ - la)[sW{l)Y{k -I)- Y{l)W{k - I)] 

w{k) = I ■ ■ ■ Ef:=J 



■E 



i„=0 



{ka - la)W,os - I) 

w{k) = -iY!i:=o---T!i:zo 

{ka-la)W,i^ {l)Y{k-l) 



i„=o 



;„=o 



For the self-referential recurrences, note that VF(0) is always determined by 
applying the original function to the zeroth-order coefficients of the respective 
power series. In addition, any terms which would cause W{k) to depend on 
itself should be omitted; their appearance is only for notational convenience. 
For some of the entries it is necessary to pick an index, a, for which the s\im 
only goes from to fca — 1. The choice of this index is arbitrary, although it 
generally must be such that fca ^ 0. 
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fix 
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Evolving fJ,f 'Kf^'^ 
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f{x,t + dt),f(x,t + dt) 
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Figure 1: This diagram illustrates the stacking technique for a second-order- 
in-time PDE using the pairwise spatial-order evolution. /^"^ indicates the yi*'^ 
spatial derivative of /. The solid arrows indicate direct (scaled) assignment. 
The dot-dashed lines indicate computation using the interpolation procedure 
detailed in Appendix |Xj The label "R" indicates computation using the recur- 
rence relation derived from the PDE. The right-curly bracket indicates the linear 
combinations from taking derivatives of the Taylor expansion which defined the 
recursion coefficients. 
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Figure 2: The relative error compared to the analytic solution for the system 
described in Section [4?2] with 1 DTM coefScient per point per temporal order. 
The interpolation uses a six-grid-point radius (13 points per neighborhood). 
CoefRcients of order 13 were the highest-order coefficients used. Also shown is 
the relative error using an interpolation neighborhood radius of seven (15 points 
per neighborhood). 
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Figure 3: The relative error compared to the analytic solution for the system 
described in Section [4?2] with 2 DTM coefficients per point per temporal order. 
This plot shows cases where the interpolation uses the two highest orders from 
all points within a three-grid-point radius (7 points per neighborhood), a four- 
grid-point radius (9 points per neighborhood), and a five-grid-point radius (11 
points per neighborhood). Coefficients of order 14 were the highest-order coeffi- 
cients used. For comparison, also shown is the case where only the zeroth-order 
components are stored as a seven-grid-point interpolation radius is used. As can 
be seen, using two orders per grid point yields a much-more-accurate evolution 
than using only one. 
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Figure 4: The relative error compared to the analytic solution for the system 
described in Section [T2| with 14 DTM coefficients per point per temporal order. 
This plot shows cases where the interpolation uses a five-grid-point radius (11 
points per neighborhood). Each pair of spatial orders is evolved together, and, 
because the PDE is linear, independently of all of the other orders. Coefficients 
of order 25 were the highest-order coefficients used (coefficients of order 14 were 
the highest-order coefficients used for the first pair, and so on). The neighbor-to- 
neighbor constraint-violation error, also shown, becomes less than the analytic 
error. 
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Figure 5: The relative error compared to the analytic solution for the system 
described in Section 4.3 with 2 DTM coefficients per point per temporal or- 
der. This plot shows cases where the interpolation uses a five-grid-point radius 
(11 points per neighborhood). Coefficients of order 17 were the highest-order 
coefficients used. The RK4 scheme shown used N — 156, dx = 0.28125 and 
dt = 0.001/2650 = 3.77358e - 07, but using dt = 0.001 with the RK4 scheme 
yields an almost-identical curve. 
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Figure 6: The relative error compared to the analytic solution for the system 
described in Section [473| with 17 DTM coefficients per point per temporal or- 
der. This plot shows cases where the interpolation uses a five-grid-point radius 
(11 points per neighborhood). Each pair of spatial orders is evolved together. 
Coefhcients of order 32 were the highest-order coefficients used (coefficients of 
order 17 were the highest-order coefficients used for the first pair, and so on). 
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Figure 7: The relative error compared to the analytic solution for the system 
described in Section 4.3 with 17 DTM coefficients per point per temporal order 
computed using the double-double type from Bailey's qd package [T^ , which has 
twice the working precision of an IEEE double. This plot shows cases where 
the interpolation uses a five-grid-point radius (11 points per neighborhood). 
Each pair of spatial orders is evolved together. Coefficients of order 32 were the 
highest-order coefficients used (coefficients of order 17 were the highest-order 
coefficients used for the first pair, and so on). 
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